The dilemma of underestimating freshwater biodiversity: morphological and molecular approaches

Background Anthropogenic impacts on freshwater habitats are causing a recent biodiversity decline far greater than that documented for most terrestrial ecosystems. However, knowledge and description of freshwater biodiversity is still limited, especially targeting all size classes to uncover the distribution of biodiversity between different trophic levels. We assessed the biodiversity of the Lower Rhine and associated water bodies in the river’s flood plain including the river’s main channel, oxbows and gravel-pit lakes, spanning from the level of protists up to the level of larger invertebrate predators and herbivores organized in size classes (nano-, micro, meio- and macrofauna). Morphological diversity was determined by morphotypes, while the molecular diversity (amplicon sequencing variants, ASVs) was assessed through eDNA samples with metabarcoding targeting the V9 region of the 18S rDNA. Results Considering all four investigated size classes, the percentage of shared taxa between both approaches eDNA (ASVs with 80–100% sequence similarity to reference sequences) and morphology (morphotypes), was always below 15% (5.4 ± 3.9%). Even with a more stringent filtering of ASVs (98–100% similarity), the overlap of taxa could only reach up to 43% (18.3 ± 12%). We observed low taxonomic resolution of reference sequences from freshwater organisms in public databases for all size classes, especially for nano-, micro-, and meiofauna, furthermore lacking metainformation if species occur in freshwater, marine or terrestrial ecosystems. Conclusions In our study, we provide a combination of morphotype detection and metabarcoding that particularly reveals the diversity in the smaller size classes and furthermore highlights the lack of genetic resources in reference databases for this diversity. Especially for protists (nano- and microfauna), a combination of molecular and morphological approaches is needed to gain the highest possible community resolution. The assessment of freshwater biodiversity needs to account for its sub-structuring in different ecological size classes and across compartments in order to reveal the ecological dimension of diversity and its distribution. Supplementary Information The online version contains supplementary material available at 10.1186/s12862-024-02261-y.


Background
As an international stream, the River Rhine is the largest federal waterway in Germany and of corresponding socio-economic importance.Here, clusters of rocks, so called ripraps, can often be found as a result of river straightening [1].Moreover, the demand of earth-bound resources, of which gravel is particularly available in the grounds of flood plains and thus mined, has furthermore resulted in the creation of artificial lake systems along the river's course.Such anthropogenically changes in the cultural landscape of Central Europe within the past decades' and our intensive utilization of freshwater habitats caused degradation, pollution, water extraction, foodweb disturbance and invasive species introduction into freshwater habitats.This is causing a recent biodiversity decline far greater than in most terrestrial ecosystems [2,3].Flowing waters and the associated limnic ecosystems of floodplains are particularly impacted by the biodiversity decline and other aspects of global change [4].Despite the imminent threat to freshwater ecosystems and their socio-ecological importance, knowledge on freshwater biodiversity is still limited.However, knowledge of the species diversity and richness is essential for generating a better understanding of freshwater ecological processes which in turn provide crucial ecosystem services fundamental for human livelihoods and wellbeing [2].
Many freshwater biodiversity studies have been conducted on one or a few specific taxonomic groups [1,5,6].Certain taxonomic groups within e.g. the macrozoobenthos as well as phytoplankton, have been integrated into biodiversity monitoring programs, where they are routinely surveyed to evaluate factors such as the overall ecological status of freshwater ecosystems [7][8][9].However, investigations that integrate across all size classes are scarce [10,11], although numerous studies in different frameworks and habitats have already demonstrated the complex interactions between organisms within the different trophic levels, as well as between size classes [10,12].While top-predators like fish are often known to regulate the population of smaller size classes, protists are shown as the essential group for the transfer of energy by connecting lower and higher trophic levels, consequently affecting the whole food web [10,12].
Given the diversity of species within the different size classes in freshwater systems, biodiversity assessments at habitat-scale and for the magnitude of different freshwater ecosystems is unrealistic to be fulfilled in projects of reasonable size and duration when classical morphology-based identification is considered exclusively.While enabling quantitative assessments of a community, morphological methods are time-consuming because sampling techniques vary in suitability depending on the different taxa and size classes [13][14][15].Furthermore, correct taxon determination based on morphological characteristics requires a notable amount of knowledge.It is necessary to include at least one expert for each size class to sufficiently cover the community's diversity.However, biodiversity monitoring is already shifting from a traditional morphological to a genetic based approach for the registration of species.Environmental DNA (eDNA) metabarcoding approaches have been identified amongst the most promising approaches to overcome previous limitations of solely morphological approaches, by their ability to detect an almost complete proportion of the species diversity in a given environmental sample [16][17][18][19][20].Not only do they have the potential to significantly decrease time, costs and required knowledge but eDNA sequencing rather detects too many than too few taxa inhabiting an area [13][14][15][18][19][20][21].However, the inability to differentiate between life stages, intact and living organisms, ingested, or extraneous tissue, to calculate abundances and biomass, or to exclude taxa through unfitted barcodes are remaining obstacles with metabarcoding [16,[22][23][24].Overall, it is particularly relevant to rely on standardized procedures for biodiversity monitoring, which are as less error prone as possible, time and cost efficient and independent from exclusive expert knowledge.Thus, a combination of multiple techniques is needed, especially when targeting all size classes at once.
The aim of this study was to investigate species diversity from protists up to larger invertebrates organized in size classes (nano-, micro, meio-and macrofauna) from a long-term ecological research site located in Germany at the Lower Rhine in North Rhine-Westphalia including two riprap river Rhine sections, two oxbows as well as two gravel-pit lakes in the floodplain.The aim of the LTER-D project REES is to investigate the eco-evolutionary dynamics along a trophic cascade, integrating species representatives with regards to molecular evolution as well as size, structure and distribution patterns of populations over time and space.Within this study the biodiversity of the targeted REES habitats was initially assessed to obtain a status quo of the water bodies and to create a biodiversity inventory that allows for the choice of representative species along trophic levels for the future development of the long-term study.Diversity was assessed through morphotype richness and metabarcoding, as well as accounting for sediment composition.Furthermore, this study gives insight into how much time, expertise and what kind of mismatch can be expected by targeting biodiversity from a morphological and molecular perspective when people of different expertise and career level are contributing to the assessment.

Sampling site
The study area lies within the LTER-D REES site [25] at the Lower Rhine in North Rhine-Westphalia (district Rees) in Germany, which includes the original flood plain area of the River Rhine (partly separated by dikes) including several gravel pit lakes, Rhine oxbows and abandoned meanders, as well as the main river.Littoral benthic samples of four standing water bodies within this area (oxbow Bienen, oxbow Grietherort, gravelpit lake Reeser Meer Norderweiterung, gravel-pit lake Reeser Meer Süd) and thick epilithic biofilm communities of two ripraps in the River Rhine (Grietherort and Cologne) were collected between 8th of July and 1st of September 2021 (Table 1, Fig. 1).
The two oxbow sites have different connections to the main river.While oxbow Bienen (site A) is an old oxbow only connected to the River Rhine at extreme flood events, oxbow Grietherort (site B) is regularly connected to the main river.Both areas are under nature conservancy, whereby oxbow Bienen (site A) is particularly distinguished as a bird sanctuary.Both sites' shores showed occasional trees and grasses, accumulations of water lilies (Nymphaea), common reed (Phragmites) as well as nearby stinging nettles (especially site B).Additionally, site A's vegetation contained further patches of reedmace (Typha), mint (Mentha), as well as duckweed (Lemna) in and around the oxbows shore.
The two riprap sites were dominated by basalt boulders ranging from approximately ten to 45 cm in diameter.The riprap in Cologne (site F) and the riprap in Grietherort (site E) were located at the impact slope of the Rhine.Besides thick biofilm layers overgrowing each submerged rock, no notable macrophytes were observed below water level.
Regarding the two groundwater-fed gravel-pit lakes, Reeser Meer Norderweiterung (site C) and Reeser Meer Süd (site D), the habitats' shores contained more soft sediment than big boulder, while site D showed several rock accumulations above and below the water level.Both lakes are still being used for gravel extraction.With several kilometers distance to the River Rhine and separation by the main dike, both sites' structural features are unaffected by flooding events; however, the connection via groundwater changes in correspondence with a changing water level of the River Rhine, too.
Both sites shared aquatic patches of pondweed (Potamogeton) as well as stoneworts (Chara).Both, together with the waterweed Elodea were especially dominant at site D's shore region.While several patches of Phragmites grew near site D's shore only few individuals were observed at the site of gravel-pit lake C. A unique characteristic of site C is the presence of a dense population of the great crested newt (Triturus cristatus) and the presence of only one fish species, the sunbleak (Leucaspius delineatus), which was first observed in 2019.Before that, the Reeser Meer Norderweiterung had been fish free.

Sampling procedure
At each sampling site, nanofauna (2-20 µm), microfauna (20-200 µm), meiofauna (200 µm-2 mm), and macrofauna (< 2 mm) [27] samples were taken in addition to eDNA samples for metabarcoding analysis.Macrofauna of the lentic water bodies was sampled in triplicates using a 25 × 25 cm benthos net (mesh size 500 µm) covering 5 m 2 (oxbow site A) or 0.53 m 2 (oxbow site B, gravel-pit lakes site C and D) per replicate.The net was pulled across the bottom, so that the organisms of the overlaying water and top 2 cm of sediment ended up in the net.The sampling area was marked up with lines and disturbance of the area prior to the sampling was avoided.For the other three size classes, sampling corer (diameter 4.4 cm) were used to collect 3.7 cm deep sediment cores together with 7 cm (3.7 cm for site A, oxbow Bienen due to the low water level) of overlaying water.Ten cores per triplicate (final volume 1000 ml) were filtered through a 500 µm net to remove large particles, leaves and macrofauna.The mixture of one replicate was evenly transferred into four 250 ml plastic beakers (subsamples).For metabarcoding studies, one subsample (250 ml) was filtered through cellulose nitrate filters with a pore size of 0.45 µm (Sartorius Stedim Biotech, Sartorius AG, Göttingen, Germany).Filters were immediately fixated with the salt solution DESS [28] in 50 ml tubes.The other three subsamples (250 ml each) were individually filtered through a meiofaunal sieve (pore size 44 µm), while collecting the filtered water for diversity analyses of protists belonging to the nano-and microfauna.The retained meiofauna on the sieves was suspended in filtered water and transferred into 50 ml tubes.One subsample was fixated with formaldehyde for qualitative analysis, one was fixed in DESS for molecular analyses, and the third was left unfixed for quantitative analyses in the laboratory.
In case of the riprap samples, the benthic community of the submerged rocks was brushed off and sucked in by a pond vacuum cleaner (PONDOVAC 4, OASE GmbH, Hörstel, Germany).To sample eDNA as well as size-classes smaller than macrofauna, biofilms were collected along a transect of defined length with a total surface area of 0.017 m 2 .The area to be sampled was marked up with the help of calibrated ropes laying on the stones.The collected suspension (~ 1 L) was processed identically as the previously described sediment samples.For the macrofauna, larger amounts of biofilm were sampled to increase the possibility to collect the entire benthic community.Organisms were retained by a 500 µm net placed over the vacuum's exit pipe and the volume of run-through water was determined.Sampling continued until approximately 30 L to 40 L were collected.

Grain size measurement
Two sediment samples (cores) of the oxbows and gravelpit lakes were taken to investigate sediment grain size composition of the upper 3.7 cm.Organic matter like leaves or branches were discarded and the wet weight of the sediment was measured.Each sample was dried overnight in a compartment dryer at 60 °C (Memmert GmbH + Co. KG, Schwabach, Germany) and the total dry weight was measured.Different sieves (63 µm, 125 µm, 250 µm, 0.5 mm, 1 mm, 2 mm) were used to determine the relative contribution of each size fraction to the total dry weight.By subtracting the weight of these size fractions from the total dry weight, the amount of fine sediment (< 63 µm) was calculated and the sediment fractions associated to the sediment types gravel (> 2 mm), very coarse (2 mm -1 mm), coarse (1 mm -0.5 mm), medium (500 µm -250 µm), fine (250 µm -125 µm) and very fine sized sand (125 µm -63 µm), as well as silt together with clay (< 63 µm) [29].

Metabarcoding
DESS preserved sediment samples were vortexed for two minutes and centrifuged (4000 × g for 20 min at 4 °C, Megafuge 2.0 R, Heraeus Instruments).Environmental DNA was then extracted from 1 g sediment of each replicate sample (a total of 3 g per site) using the DNeasy PowerLyzer PowerSoil DNA isolation kit (Qiagen, Hilden, Germany) according to the manufacturer's protocol.Prior to the application of the kit, sediment samples were washed with three washing solutions to improve the success of DNA amplification by PCR [54].Total DNA was quantified using Quantus ™ Fluorometer (Promega, Wisconsin, USA).PCR amplifications of the hypervariable V9 region of the 18S rDNA gene was performed with 12.5 µl of the 2X VWR Red Taq DNA polymerase Master Mix (VWR, Erlangen, Germany) and the forward/reverse primer-pair 1389F (5′-TTG TAC ACA CCG CCC-3′) and 1510R (5′-CCT TCY GCA GGT TCA CCT AC-3′) [55].The PCR mixtures (25 µl final volume) contained 30 ng of total DNA template per site (10 ng per replicate) with 0.35 µM final concentration of each primer.PCR amplifications (98 °C for 30 s; 25 cycles of 10 s at 98 °C, 30 s at 57 °C, 30 s at 72 °C; and 72 °C for 10 min) of all samples were carried out with a reduced number of cycles to avoid the formation of chimeras during the plateau phase of the reaction.Amplifications per site were conducted in triplicates in order to smooth the intra-sample variance while obtaining sufficient amounts of amplicons for Illumina sequencing.PCR products were checked on a 2% agarose gel for amplicon lengths and fragment sizes were determined by comparison with a 50 bp DNA Ladder (Nippon Genetics Europe, Düren, Germany).Amplicons were then pooled and purified using the FastGene Gel/ PCR Purification Kit (Nippon Genetics Europe, Düren, Germany).For subsequent quality measures during data analysis for protists, we used an in vitro community ("mock community") comprising DNA of nine different protist cultures belonging to representatives of the main protist supergroups [56].We used this mock community to correct for the occurrence of non-reliable sequences.We determined the specific minimum threshold for the samples based on the corresponding mock community.This threshold was calculated by considering the proportion of the lowest read count of an ASV in the mock community dataset that could be attributed to the cultured species.ASVs in the sample dataset with read counts lower than this calculated proportion (0.057%) were excluded.Similar thresholds were calculated in previous studies using the same mock community [56,57].Bridge amplification and paired-end (2 × 150 bp) sequencing of the amplified fragments were performed using a NovaSeq 6000 system at the Cologne Center of Genomics (CCG) in Germany.

Raw sequencing data processing and taxonomic assignment
After sequencing, the raw reads were demultiplexed (parameters: -no-indels, -overlap 8, e 0) and barcode and primer sequences (parameters: -no-indels, -m 30, e 0.2) were clipped using cutadapt version 2.8 [58].Amplicon sequencing variants (ASVs) were generated by using DADA2 version 1.26 [59] implemented in R v.4.3.0 [60].Filter and trimming parameters were set to maxEE = 1, truncQ = 11, truncLen = (125, 120), and maxN = 0 for quality filtering of the reads, followed by training of error models using the learnErrors function of DADA2 with default settings.Dereplication of sequences was done with the derepFastq function and the inference of ASVs was done with the dada function and default settings.Paired reads were merged with a minimum overlap of 12 nucleotides with mergePairs.Chimeric sequences were removed using the removeChimeraDenovo function.We used the PR 2 database version 4.14.0 [61] for taxonomic assignment of ASVs via the pairwise alignment function usearch_global of VSEARCH v2.18.0 [62].All unassigned sequences were removed, keeping only ASVs with a pairwise identity of ≥ 80% to a reference sequence.First, we divided our dataset into a protist and metazoan dataset.Since this study focuses on zooplankton/zoobenthos, we excluded phototrophic protist sequences (determined on the basis of taxonomic assignment) for the protist dataset including Ochrophyta, Haptophyta, Filosa-Chlorarachnea within the Cercozoa as well as the Cryptomonadales and Pyrenomonadales within the Cryptophyta.Moreover, we removed Metazoa, Archaeplastida and Fungi for the protist dataset.As described by Dünn and Arndt [56], as a last filtering step, we used a minimum threshold of the previously described mock community by calculating the proportion of the lowest read number of an ASV in the mock community data set that could be assigned to the cultured species with a pairwise identity of 100%.ASVs in the sample data sets with a smaller read number than this calculated proportion were discarded.For the metazoan dataset, besides ASVs assigned to Craniata, we removed terrestrial and marine groups from the final ASV dataset as these might be the result of database limitations.The protist dataset was divided into nanoand microfauna and the metazoan dataset into meioand macrofauna based on size ranges of taxa found in literature.

Community diversity analysis
Statistical analyses were conducted with R v.4.3.0 [60].Binary-Jaccard distances were used as a measure of beta-diversity by calling the function vegdist within the "vegan" package [63].The Jaccard distance values were then used for the unweighted pair-group method with arithmetic means (UPGMA) cluster analyses (hclust function).Results of the cluster analyses were visualized in dendrograms by using "ggplot2" [64].Non-metric multidimensional scaling (NMDS) was performed with Jaccard distances using "vegan" to show variation of habitats based on morphotype detection.Observed mean morphotype richness per size classes of respective sampling sites were compared amongst each other via Kruskal-Wallis one-way analysis of variance and Dunn's post hoc analysis with three replicates per size-class and site (n = 72).Venn diagrams [65,66] were used to visualize the number of shared and unique ASVs and morphotypes between the three habitats and the different size classes.Heatmaps were created by using the package "ampvis2" [67].For the heatmaps, only ASVs with a 98-100% sequence similarity to deposited sequences in the PR 2 database [61] were kept and clustered to a predicted genus level.With regards to the relative abundance of this ASV sequence similarity range, the 35 most abundant genera per size class were shown in the heatmaps.

Results
While investigating different size classes and taxonomic groups within this study, we looked at different levels of molecular diversity, including the molecular diversity based on ASVs with 80-100% sequence similarity to reference sequences (ASV 80-100% ) and at the diversity on a predicted genus level based on ASVs with 98-100% sequence similarity to reference sequences followed by a taxonomic clustering based on genus level (ASV 98-100% ).These ASVs 98-100% were used for a direct comparison between the morphological diversity based on the morphological approach (morphotypes) and molecular diversity on a predicted genus level based on metabarcoding.

Samples and sediments
In total, 12 samples were taken from the six different sampling sites including four size classes.Three replicates were sampled per size class.All samples could successfully be investigated according to the full spectrum of analysis of this study, except for metabarcoding of the oxbow Grietherort (site B) due to PCR failure.Alongside with the biological samples, sediments of each sampling site were analyzed to better define the habitat characteristics on the micro scale via grain size categories.
Four sites (oxbows and gravel-pit lakes) were characterized by fine sediment, oxbow site A and gravel-pit lake site D by gravelly muddy sand, the gravel-pit site C by gravel sand and the oxbow site B as muddy sandy gravel (Additional file 1, Fig. S1).

Total number of morphotypes and ASVs
This study investigates biodiversity from three different perspectives: morphological diversity (determined by morphotypes), molecular diversity (determined by ASVs with 80-100% sequence similarity to reference sequences) and the diversity on a predicted genus level (determined by ASVs with 98-100% sequence similarity to reference sequences).
Overall, zoobenthos richness summarized across all size classes and stations showed that out of the total of 191 morphotypes (determined to a lower taxon than the supergroup level) the majority belonged to macrofauna (103 morphotypes) and the lowest number of morphotypes was observed within meiofauna (22 morphotypes).For the protist communities, 24 nanofauna and 42 microfauna morphotypes were detected (Additional file 2, Table S1).

Habitat specific zoobenthos richness
The NMDS analysis based on morphotype richness showed that all three different habitats clustered separately with a slight overlap between oxbows and gravelpit lakes (Fig. 2A).The mean richness considering all size classes within each of the three habitat types was almost the same in the gravel-pit lakes and the River Rhine (42.3 ± 11.5 ind.and 40 ± 4.05 ind., respectively), while the oxbows showed the lowest mean morphotype richness with 32.5 ± 8.6 individuals (Fig. 2A Violin plot).Overall, zoobenthos richness with regard to morphotypes across all size classes was not significantly different between habitats (Kruskal-Wallis chi-squared = 3.03, df = 2, p-value = 0.22).
The zoobenthic community of the River Rhine was different from all other habitats with 64 unique morphotypes, explaining 34.4% of the community richness (Fig. 2B).Oxbows inhabited the lowest number of 69 morphotypes, whereas 93 and 86 morphotypes were found in the gravel-pit lakes and the River Rhine, respectively.Only a small proportion (10 morphotypes, 5.4%) was shared between all habitat types, whereas gravelpit lakes and oxbows shared 30 morphotypes (16.1%) as indicated in the NMDS plot (Fig. 2A).

Taxonomic comparison of richness assessed via ASVs and morphotypes
Hierarchical clustering of ASVs 80-100% and morphotype richness showed that each habitat formed a separate cluster within the nano-, micro-and macrofauna (Figs.3A, C, 4A and 5A).For the meiofauna, however, the clustering differed in branching of the gravel-pit lakes, while the cluster for the River Rhine sites was consistent (Fig. 4A).Overall, nano-and microfauna richness were less similar between sampling sites (Jaccard distances between 0.56-0.93)when compared to the meio-and macrofauna (Jaccard distances between 0.39-0.63)(Figs.3A, C, 4A  and 5A).

Nanofauna
When considering all ASVs with a sequence similarity of 80-100% for the nanofauna community, ASVs 80-100% could be assigned to 32 orders.However, only morphotypes belonging to five out of these 32 classes belonging to five supergroups and Incertae sedis could be identified, including morphotypes belonging to the Filosa Sarcomonadea (Rhizaria), Cryptophyceae (Cryptista), Kinetoplastea and Euglenida ("Excavata"), Ancyromonadida (Incertae sedis) (Fig. 3A).Except for the Ancyromonadida, we also detected a high ASV 80-100% richness for these division levels, where morphotypes were found.The highest morphotype richness could be detected for Kinetoplastea and Euglenida, while especially samples from the oxbow revealed a high euglenid ASV 80-100% richness (Fig. 3A).Out of the 24 nanofauna morphotypes, only eight could be determined down to genus level, including kinetoplastids (Neobodo, Rhynchomonas), cercomonads (cf.Cercomonas), euglenids (Petalomonas, Entosiphon, Peranema), cryptophyceans (Goniomonas), ancyromonadids (Ancyromonas = Nutomonas) (Additional file 2, Table S1).When considering ASVs with a sequence similarity of 98-100% with a taxonomic clustering on predicted genus level, only three genera out of the 35 most abundant ASVs 98-100% could also be detected morphologically, however, with great discrepancies when looking at site level matches (e.g.Rhynchomonas and cf.Cercomonas) (Fig. 3B).The most abundant nanofauna genus, the cercozoan Rhogostoma, based on read abundance, could not be detected morphologically in the samples (Fig. 3B).The second most abundant genus, the kinetoplastid Neobodo, was detected in all habitats and sites by metabarcoding, but morphotype detection revealed this genus to be present only in the River Rhine and one gravel-pit lake (site D).

Microfauna
When considering all ASVs with a sequence similarity of 80-100% for the microfauna community, ASVs 80-100% could be assigned to 32 classes belonging to five supergroups.However, only morphotypes belonging to ten out of these 32 classes could be identified, including morphotypes belonging to the Ciliophora (Litostomatea, Nassophorea, Oligohymenophorea, Phyllopharyngea, Prostomatea 1, Spirotrichea), Dinoflagellata (Dinophyceae) and Amoebozoa (Fig. 3C).Ciliophora had the highest ASV 80-100% richness at all sites, with the majority of morphotypes assigned to this group (37 out of 42 total morphotypes).Within the Ciliophora except for the Nassophorea and Prostomatea 1, we also detected a high ASV 80-100% richness for these division levels, where morphotypes were found.The highest morphotype richness within the microfauna community could be detected for the Spirotrichea, while also the highest ASV 80-100% richness was observed for this group (Fig. 3A).Only four amoebozoan morphotypes were found, although the number of ASV 80-100% and richness was much higher.While we recovered ASV 80-100% for several divisions within the rhizarians, we were not able to detect them morphologically.Out of the 42 microfauna morphotypes, 21 could be determined down to genus/species level, including species belonging to the Spirotrichea, Phyllopharyngea, and Oligohymenophorea (Additional file 2, Table S1).When considering ASVs with a sequence similarity of 98-100% with a taxonomic clustering on In several cases reference sequences were not assigned to genus level, thus, a higher taxonomic level is shown.Numbers within the heatmap correspond to the number of ASVs assigned to this genus per site.Taxonomic groups correspond to class and genus levels in the PR 2 database classification.The sequential color code corresponds to the relative abundance of reads with a sequence similarity of 98-100% assigned to the respective genus to either the nano-or microfauna size class.Pink asterisks indicate, if genera could also be detected morphologically.Nanofauna protist drawings are adapted from literature [31,68].Microfauna protist silhouettes are from PhyloPic [69] contributed by Guillaume Dera, 2023 (CC0 1.0 Universal Public Domain Dedication) and Yan Wong, 2013 (Attribution 3.0 Unported, https:// creat iveco mmons.org/ licen ses/ by/3.0/) predicted genus level, only two genera out of the 35 most abundant ASVs 98-100% could also be detected morphologically, namely the genera Trithigmostoma and Euplotes (Fig. 3D).

Comparison of methods
The performance and outcomes of biodiversity surveys have been examined within several studies with a particular focus on the comparative effectiveness of morphology and DNA-based approaches of specific taxonomic groups highlighting discrepancies that may arise in the application of these methods [70][71][72][73][74][75][76].In addition, single primer pairs or combinations of primer pairs sequencing the same or different regions have been used within metabarcoding studies targeting different taxa [77][78][79][80][81].The choice of a primer pair is dependent on the research question, the targeted taxonomic groups and the sequencing technique (e.g.short vs. long reads).For protists the 18S rRNA is used as a universal marker to assess protist biodiversity targeting mainly the hypervariable V4 or V9 regions (~ 420 bp and ~ 130 bp, respectively) in metabarcoding studies [55,78,82,83].For several protist taxa the V9 region can be used to detect down to the species level [84], while for others this region was only sufficient to distinguish down to the genus level [85].A In several cases reference sequences were not assigned to genus/species level, thus, the deepest taxonomic level such as the order is shown.Numbers within the heatmap correspond to the number of ASVs 98-100% assigned to these genera for each of the five sites.Taxonomic genus rank corresponds to genus levels in the PR2 database classification.The sequential color code corresponds to the relative abundance of reads assigned to the genus and is relative to all reads assigned to the macrofauna size classes.Pink asterisks behind the taxonomic names indicate, if morphotypes down to genus level could also be identified.Pink asterisks indicate, if genera could also be detected morphologically.Organism silhouettes are from PhyloPic [69].Contributed by Nathan Jay Baker, 2022 (Branchiopoda), Kamil S. Jaron, 2022 (Insecta), Tauana Cunha, 2021 (Gastropoda) and Guillaume Dera, 2023 (Rhabditophora, Trematoda, Porifera) under license CC0 1.0 Universal Public Domain Dedication.Contributed by Denis Lafage, 2020 (Malacostraca), Katie Collins, 2020 (Bivalvia) and Mathilde Cordellier, 2020 (Arachnida) under License Attribution 3.0 Unported (https:// creat iveco mmons.org/ licen ses/ by/3.0/).Contributed by B. Duygu Özpolat, 2016 (Hirudinea, Oligochaeta) under license Attribution-NonCommercial-ShareAlike 3.0 Unported (https:// creat iveco mmons.org/ licen ses/ by-nc-sa/3.0/) reliable taxonomic assignment of nematode species is based on the 18S rRNA together with the highly conserved mitochondrial region cytochrome oxidase COI [86].The hypervariable V9 region has also been used to investigate the overall metazoan community composition (e.g.[87]).While for a variety of metazoan species the COI is used as a reliable marker gene for identification, several difficulties challenging distinct species identification using COI as the marker gene were addressed by the barcoding community in the past [88][89][90], such as the "barcoding gap" [91].While targeting four size classes, we used the 18S V9 region as marker region to cover a broad spectrum of organisms of all size classes.
In several cases the precision and taxonomic resolution of biodiversity surveys had been increased by eDNA surveys [73,92,93].When considering four different size classes, the percentage of shared taxa between both approaches eDNA (ASVs 80-100% ) vs. morphology (morphotypes) was always below 15% (5.4 ± 3.9%) within our study.A study targeting phytoplankton, zooplankton and macroinvertebrates also showed a low number of taxa shared between both approaches with 7-9% [11].When we considered only ASVs with a sequence similarity of 98-100%, the number of shared taxa reached up to 43% (18.3 ± 12%) (see Additional file 3, Table S4).Different taxa were most likely recorded under the same morphotype, likely underestimating the true diversity of each site.In addition, DNA has the ability to persist for different time spans in the environment even after organisms have disappeared.In lentic water bodies this can also lead to the downstream transport of eDNA molecules beyond the natural distribution range of the living species.Extracellular DNA can be preserved from hours to days in the water column and from decades to centuries in the sediment [94].This DNA persistence can lead to false positive results that indicate the presence of organisms that are no longer part of the ecosystem being studied.The River Rhine spans several hundred kilometers of length and thus provides the highest potential for persisting eDNA molecules of upstream occurring species that might then be detected in the lower reaches, such as the sampling sites of our study.However, the molecular diversity of the Rhine samples did not outcompete the other water bodies in neither of the size classes.The overestimation of the diversity of the River Rhine on the level of ASVs due to a high influx of persistent eDNA molecules from upstream is thus rather unlikely.
Overall, many major taxonomic groups in all size classes could not be recovered with the morphological approach and on the other hand, many ASVs 98-100% of all size classes have not been found within our morphological detection (see Figs. 3B, 4B and 5B).But there were also examples of matches of both applied methods, genus-clustered ASVs 98-100% and morphotypes, that were observed within all size classes, e.g. for the macrofauna the water louse Asellus, the mysid Limnomysis as well as the dragonfly Libellula, and for the meiofauna the gastrotrich Chaetonotus.Especially for protists (nano-and microfauna), a combination of molecular and morphological approaches is needed to gain the highest possible community resolution.Within the microfauna, for example, the phyllopharyngean cyrtophorid ciliate genus Trithigmostoma, a common freshwater ciliate, was detected within both ripraps by both approaches.High read abundances of the cercozoan species Rhogostoma minus were detected by metabarcoding within the two riprap sites (E and F), but was also detected in the other two gravel-pit lakes as well as the oxbow.This delicate thecofilosean, however, could not be identified with the morphological approach during the sampling period, but was identified in the Rhine in earlier studies (H.Arndt, unpubl.).Another example is the ciliate genus Stentor, which is common in the River Rhine at Cologne [5] and was highly abundant with regards to reads from the metabarcoding approach, but could not be detected by our morphological approach at the sampling date.
One advantage of metabarcoding is that it is less reliant on taxonomic expertise [16,71] including many taxa currently not identifiable by expert taxonomists [95][96][97].However, metabarcoding can only be as precise as its database [18,22,24].Not only must a database include enough correctly assigned taxa to potentially find new invasive species, the aspect that many public databases are error-infected makes the availability of well-curated databases an important prerequisite [24].Several morphotypes identified to species level in this study have been missing in the database.Within the microfauna, one example is the detected ciliate genus Aspidisca for which sequences could be found in the used reference database (e.g. A. steini, A. aculeata, A. magna), but not for the common species identified by the morphological approach (A.cicada, A. turrita, A. lynceus).
For the macrofauna, we found Dikerogammarus morphotypes only to be present in the samples from the Ripraps in the River Rhine, which are seen as an emerging ecosystem with new substrate-specific combinations and can be dominated by non-native species such as Dikerogammarus villosus/haemobaphes [1,98,99].However, metabarcoding only indicated the presence of Gammaridae in the other two habitats (oxbows and gravel-pit lakes), namely Gammarus pulex (ASV with 100% sequence similarity) and Gammarus tigrinus (ASV with 100% sequence similarity).While no V9 sequence for Dikerogammarus was present in the used reference database, reblasting these two Gammarus ASVs (read abundance of 394 and 55) at GenBank showed a 100% similarity to many species belonging to the family of Gammaridae, but only a 94.46% similarity to Dikerogammarus species.

Habitat specific communities and biodiversity inventory
Irrespective of the size class, we observed habitat specific zoobenthos communities in each of the three investigated habitat types, with few taxa or taxonomic groups overlapping.The gravel-pit lakes showed the highest overall molecular diversity (ASV 80-100% ) of zoobenthos, while the overall morphological diversity was similar for the gravel-pit lakes and the River Rhine.We found the highest unique morphotype richness in the ripraps of the River Rhine with a distinct community composition when compared to the other habitats.While meiofauna morphotypes, especially nematodes, were underrepresented due to the needed taxonomic expertise to identify genera or species morphologically, we recovered a high genetic diversity of nematodes belonging to the order Monhysterida and Chromadorea as well as Chaetonotida with our metabarcoding approach (ASV 80-100% ) within the gravel-pit lakes.Only for the macrofauna, the morphotype richness was highest in the gravel-pit lakes, while for the other size classes (i.e. the nano-and microfauna) the highest unique morphotype richness was documented from the River Rhine.
Differences in sediment composition and macrophytes at each site might strengthen the differences in community richness between those two gravel-pit lakes.Moreover, only the sunbleak Leucaspius delineatus was present that at the Reeser Meer Norderweiterung (site C), while the Reeser Meer Süd (site D) is inhabited by many fish species including the European perch (Perca fluviatilis), the common roach (Rutilus rutilus) and the northern pike (Esox lucius).While patches of macrophytes in the shore area of lakes are generally known to provide distinct ecological niches and food sources benefitting the zoobenthos diversity [100][101][102], the high occurrence of macrophytes (Potamogeton, Chara, Elodea), especially at gravel-pit lake site C's (Reeser Meer Norderweiterung) shoreline could be one explanation for the higher ASV 80-100% within the gravel-pit lakes when compared to the River Rhine.
With this study, we tried to obtain a status quo of the water bodies that distinguishes biodiversity across the ecological size classes.With this, we display the distribution of diversity across ecological dimensions with nano-and meiofauna explaining a major proportion of the diversity that can be best assessed via metabarcoding.This ecological partitioning of diversity, however, is in stark contrast to the availability of reference information in sequence databases and highlights the need of genetic assessments of freshwater biodiversity to expand reference libraries.The biodiversity inventory of the here presented different water bodies, furthermore builds the necessary prerequisite to identify representative species along trophic cascades to aid future studies that investigate functional aspects of the freshwater communities.

Perspectives for biodiversity monitoring
Biodiversity assessments based on metabarcoding might not replace the traditional morphological approach until sequence availability in databases has been extended dramatically, but is suggested to be used as a complementary tool for biodiversity monitoring [72,73,103,104].The combination of morphotype detection and metabarcoding was particularly relevant for the biodiversity assessment of bodonids present in all three investigated habitats (genera Neobodo, Rhynchomonas, Dimastigella).This order of kinetoplastids is known to harbor a very high degree of genetic diversity as compared to whole orders of higher eukaryotes suggesting large numbers of cryptic individual species within bodonids [105].
While traditional morphological approaches as well as metabarcoding both have their advantages and limitations, metagenomics might promise a new approach for biodiversity assessment to overcome the limitations of metabarcoding such as primer dependency and PCR amplification biases.The present and further acquisition of reference genomes will help discover the genetic mechanisms underlying organisms' responses to their natural environment and with this support conservation efforts [106].Reference genomes from a wide range of species are required to map genomic variability and ultimately contribute to the conservation of genetic diversity.Various international initiatives aim to generate reference genomes representing global diversity [107].With high quality reference genomes, genomic diversity of individuals from the same species can be cost-efficiently unraveled by resequencing and aligning against it [108].

Fig. 2
Fig. 2 Richness across habitat types.A Non-metric multidimensional scaling (NMDS) analysis of the Jaccard distance matrix computed from morphotypes (n = 191 morphotypes) at the six sites (n = 3 replicates), associated to three different habitat types (Oxbow, Gravel-pit lakes, River Rhine).The top inset displays the distribution of morphotype richness for the three water body communities.The black dots and bars within violin plots represent means and SDs (replicates per habitat type n = 6).B-C Venn diagrams showing the number of unique and shared (B) morphotypes and (C) ASVs 80-100% between the three different habitat types including all four size classes (nano-, micro-, meio-and macrofauna)

(Fig. 3
See figure on next page.)Distribution and community composition of freshwater nanofauna (A-B) and microfauna (C-D).A and C Dendrogram cluster showing the similarity (Jaccard index) of size class communities of the five sediment samples in regard to species richness based on incidence-based data (presence/absence) using UPGMA clustering.Number of freshwater ASV 80-100% (grey bars) and number of morphotypes (pink bars) related to the major taxonomic protist groups are shown.Taxonomic groups correspond to class level in the PR 2 database classification.B and D Heatmap of nano-and microfauna ASVs 98-100% with a 98-100% sequence similarity and taxonomic clustering on predicted genus level.Shown are the first 35 most abundant genera (out of 43 nanofauna genera and 79 microfauna genera) with class and genus level.

Fig. 4
Fig. 4 Distributional patterns and community composition of freshwater meiofauna.A Dendrogram cluster showing the similarity (Jaccard index) of heterotrophic protist communities of the five sediment samples in regard to species richness based on incidence-based data (presence/ absence) using UPGMA clustering.Number of freshwater ASVs 80-100% (grey bars) and number of morphotypes (pink bars) related to the major taxonomic groups are shown.B Heatmap of meiofauna ASVs with a 98-100% sequence similarity and taxonomic clustering on predicted genus level.Shown are the first 35 most abundant predicted genera (out of 67 genera) with class and genus level.In several cases reference sequences were not assigned to genus/species level, thus, the deepest taxonomic level such as the order is shown.Numbers within the heatmap correspond to the number of ASVs 98-100% assigned to these genera for each of the five sites.Taxonomic genus rank corresponds to genus levels in the PR 2 database classification.The sequential color code corresponds to the relative abundance of reads assigned to the genus and is relative to all reads assigned to the meiofauna size classes with a sequence similarity of 98-100%.Pink asterisks indicate, if genera could also be detected morphologically.Organism silhouettes are from PhyloPic [69].Contributed by Mathilde Cordellier, 2020 (Arachnida), Birgit Lang, 2015 (Collembola), Maxime Dahirel, 2018 (Ostracoda), Michelle Site, 2014 (Nematoda) under License Attribution 3.0 Unported (https:// creat iveco mmons.org/ licen ses/ by/3.0/).Contributed by Siel Wellens, 2019 (Copepoda), T. Michael Keesey, 2013 (Branchiopoda), Scott Hartmann, 2013 (Gastrotricha) and Levi Simons, 2023 (Cestoda and Rotifera) under license CC0 1.0 Universal Public Domain Dedication.Contributed by Julie Blommaert, 2020 (Rotatoria) under license Attribution-ShareAlike 3.0 Unported.Contributed by B. Duygu Özpolat, 2016 (Oligochaeta) under license Attribution-NonCommercial-ShareAlike 3.0 Unported

(Fig. 5
See figure on next page.)Distributional patterns and community composition of freshwater macrofauna.A Dendrogram cluster showing the similarity (Jaccard index) of macroinvertebrate communities of the five sediment samples in regard to species richness based on incidence-based data (presence/ absence) using UPGMA clustering.Number of freshwater ASVs 80-100% (grey bars) and number of morphotypes (pink bars) related to the major taxonomic macroinvertebrate groups are shown.B Heatmap of macrofauna ASVs with a 98-100% sequence similarity and taxonomic clustering on predicted genus level.Shown are the first 35 most abundant genera (out of 115 genera) with class and genus level.

Table 1
List of the six sampling sites including sampling dates and coordinates